library(nloptr)
library(corpcor)

exp_mov <- read.csv("daily_expected_mov.csv",na.strings=c(".", "NA", "", "?"), strip.white=TRUE, encoding="UTF-8")
exp_vol <- read.csv("daily_expected_vol.csv",na.strings=c(".", "NA", "", "?"), strip.white=TRUE, encoding="UTF-8")
hist_ret<- read.csv("daily_hist_ret.csv",na.strings=c(".", "NA", "", "?"), strip.white=TRUE, encoding="UTF-8")

numrow <-nrow(exp_mov)
numcol <-ncol(exp_mov)-1

exp_mov<-exp_mov[,2:ncol(exp_mov)]
exp_vol<-exp_vol[,2:ncol(exp_vol)]
hist_ret<-hist_ret[,2:ncol(hist_ret)]
#Active Position Matrix
active_pos <- matrix()
length(active_pos) <- numrow*numcol
dim(active_pos) <- c(numrow,numcol)

eval_f <- function( x,alpha,covariance,xprime) {
  
  return( list( "objective" = -t(x)%*%alpha+0.5*t(x) %*% covariance %*% x,
                "gradient" = alpha-t(x)%*%covariance))   
}


# equalities
eval_g_eq <- function( x,alpha,covariance,xprime) {

    constr <-  rbind( sum(x))
    grad <- rbind(rep(1,numcol))
    return( list( "constraints"=constr, "jacobian"=grad ) )
}

# inequalities
eval_g_ineq <- function( x,alpha,covariance,xprime) {
  
  constr <- rbind(t(x) %*% covariance %*% x-0.1)
  grad <- rbind(2*t(x)%*%covariance)
  return( list( "constraints"=constr, "jacobian"=grad ) )
}


# initial values
x0 <- rep(0.1,numcol)
xprime=rep(0,numcol)
# lower and upper bounds of control
lb <- rep(-0.2,numcol)
ub <- rep(0.2,numcol)
local_opts <- list( "algorithm" = "NLOPT_LD_SLSQP",
                    "xtol_rel" = 1.0e-7 )

opts <- list( "algorithm" = "NLOPT_LD_SLSQP",
              "xtol_rel" = 1.0e-7,
              "maxeval" = 1000000,
              "local_opts" = local_opts )


for ( i in 1592:numrow)
{
   
    histret <- hist_ret[(i-1):(i-512),1:numcol]
    hist_cor = cor.shrink(histret)
    vol=as.numeric(exp_vol[i,])
    exp_cov=diag(vol)%*%hist_cor%*%diag(vol)
    alpha=as.numeric(exp_mov[i,])
  
    res <- nloptr( x0=x0,
                   eval_f=eval_f,
                   lb=lb,
                   ub=ub,
                   eval_g_eq=eval_g_eq,
                   eval_g_ineq=eval_g_ineq,
                   opts=opts,
                   alpha=alpha,
                   covariance=exp_cov,
                   xprime=xprime)
    xprime=res$solution
    active_pos[i,]=res$solution
}  

write.csv(active_pos,"weights.csv")